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ABSTRACT 



Aims. In the context of the core instability model, we present calculations of in situ giant planet formation. The oligarchic growth 
regime of solid protoplanets is the model adopted for the growth of the core. This growth regime for the core has not been considered 
before in full evolutionary calculations of this kind. 

Methods. The full differential equations of giant planet formation were numerically solved with an adaptation of a Henyey-type code. 
The planetesimals accretion rate was coupled in a self-consistent way to the envelope's evolution. 

Results. We performed several simulations for the formation of a Jupiter-like object by assuming various surface densities for the 
protoplanetary disc and two different sizes for the accreted planetesimals. We first focus our study on the atmospheric gas drag that 
the incoming planetesimals suffer. We find that this effect gives rise to a major enhancement on the effective capture radius of the 
protoplanet, thus leading to an average timescale reduction of ~ 30% - 55% and ultimately to an increase by a factor of 2 of the final 
mass of solids accreted as compared to the situation in which drag effects are neglected. In addition, we also examine the importance 
of the size of accreted planetesimals on the whole formation process. With regard to this second point, we find that for a swarm of 
planetesimals having a radius of 10 km, the formation time is a factor 2 to 3 shorter than that of planetesimals of 100 km, the factor 
depending on the surface density of the nebula. Moreover, planetesimal size does not seem to have a significant impact on the final 
mass of the core. 

Key words. Planets and satellites: formation - Solar System: formation - Methods: numerical 



1. Introduction 

Although the giant planets of the Solar System are known from 
centuries ago and the existence of other planetary systems har- 
bouring Jupiter-like objects was confirmed in 1995 (Mayor & 
Queloz 1995), the mechanism of giant planet formation is still 
not solved. In fact, two major models currently coexist: the core 
accretion model — alternatively referred to as the nucleated in- 
stability model — and the gas instability model (for a detailed 
review on the theory of giant planet formation see, e.g., Wuchterl 
et al. [20001 Lissauer & Stevenson |2007l Durisen et al. [2007l >. 

The core accretion model assumes that giant planets form, 
roughly speaking, in a two-step process (e.g. Mizuno 1980; 
Bodenheimer & Pollack |19861 l. First, analogous to terrestrial 
planets, a solid core is formed by coagulation of planetesimals. 
Once the core is massive enough to capture a significant amount 
of gas from the surrounding nebula, the characteristic exten- 
sive envelope of these objects forms. Gravitational gas binding 
occurs mainly in two different regimes (Pollack et al. [1996). 
Initially, gas accretion proceeds relatively quietly. This stage 
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ends when the mass of the envelope is comparable to the mass 
of the core. When this critical mass is reached, a gaseous run- 
away growth completes the formation of the planet, accreting 
large amounts of gas in a very short time. 

Evidence that the giant planets of the Solar System have a 
core approximately a dozen times the mass of Earth, favours 
the nucleated instability scenario over the gas instability hy- 
pothesis. However, one of the problems that still remains open 
with the core accretion model is how to reconcile the forma- 
tion timescales evolving from numerical simulations with those 
obtained by observations of circumstellar discs (e.g., Haisch et 
al. 1200 11 Chen & Kamp 2004). Observational evidence restricts 
the lifetime of protoplanetary discs to less than 10 7 yr. This 
fact puts a limit on the whole process: giant planets should be 
formed before the dissipation of the disc. After the pioneer- 
ing work of Pollack et al. ( 119961 ), in which detailed calcula- 
tions of giant planet formation were performed, efforts have been 
made to solve this problem. In particular, different scenarios that 
could alleviate the timescale issue have been explored, includ- 
ing: planet migration (Alibert et al. 120051 ), grain opacity reduc- 
tion (Podolak 2003 ; Hubickyj et al. [2005 ), and core formation 
in the centre of an anticyclonic vortex (Klahr & Bodenheimer 
2006) among others. To this end, most authors who adopt a 
time-dependent solid accretion rate for the core (Pollack et al. 
[1996; Hubickyj et al. 120031 Alibert et al. 120051 ) prescribe a rapid 
one (Greenzweig & Lissauer 1992) which leads to the formation 
of a massive core in a few hundred thousand years. This accre- 
tion rate slows down when most of the material within the feed- 
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ing zone is swallowed by the embryo. With this core accretion 
rate, the formation timescale of a core of about 10 M ffl usually 
only represents a small fraction of the whole formation process 
(Pollack et al. 1996). However, N-body simulations show that a 
solid embryo with the Moon mass, or even smaller, could gravi- 
tationally perturb the swarm of planetesimals around it, heating 
the disc and decreasing solid accretion long before the isolation 
mass is achieved (Ida & Makino 11993I ). The runaway growth 
of the large planetesimals (the first stage in protoplanet build- 
ing) then switches to a slower, self-limiting mode. This regime 
is known as "oligarchic growth" (Kokubo & Ida 119981 2000, 
2002; Thommes et al. 120031) . Numerical and semi-analytical cal- 
culations were made to estimate the formation timescale of a 
solid protoplanet via the oligarchic growth (Ida & Makino 1993 1 
Thommes et al. 120031 1. However, no full evolutionary calcula- 
tions of gaseous giant planet formation prescribing this accretion 
rate have been made to date. 

In the frame of the core instability model, this paper is in- 
tended to study the in situ formation of a protoplanet in a cir- 
cular orbit around the Sun, at the current position of Jupiter. 
To this end, the numerical code developed by Benvenuto & 
Brunini (2005 ) for self-consistent calculation of giant planet for- 
mation and evolution has been updated. Motivated by the work 
of Thommes et al. (2003 ), where the oligarchic regime was ap- 
plied to growing solid embryos in the absence of an atmosphere, 
we generalise this accretion model in order to adopt it as the 
time-dependent accretion rate for the core. The accretion rate 
employed in Thommes et al. ( 120031 1 is modified to take into 
account the gas drag that incoming planetesimals suffer due to 
the presence of the atmosphere, as it enlarges the protoplanet's 
capture cross-section. Inaba & Ikoma (2003 ) derived a semi- 
analytical core accretion rate for a protoplanet surrounded by a 
gaseous atmosphere and estimated the enhancement in the col- 
lisional cross-section it produced. In their study, the structure 
of the atmosphere is calculated only for certain core masses and 
not as a result of an evolutionary sequence of models. In our 
study, the solids accretion rate is coupled to the gas accretion 
rate, leading to a self-consistent evolutionary calculation of the 
envelope's structure. Hence, a more accurate estimation of the 
enlargement of the effective accretion cross-section of the pro- 
toplanet can be made. 

This paper is organised as follows: In Sect. [2] we describe 
the improvements in our numerical code, placing particular em- 
phasis on the treatment of the mass accretion rate of planetesi- 
mals onto the core. In connection with the growth of the core, 
we summarise some conceptual aspects and the main results re- 
garding the oligarchic growth of a solid protoplanet. In Sect. [3] 
we present a simulation which is compared to results obtained 
by other authors which use a different planetesimals accretion 
rate. This is intended to show the impact that the selected core 
accretion model has on the formation of a giant planet. We then 
present the results arising from running our code for simulations 
of the formation of a Jupiter-like planet for several protoplan- 
etary disc surface densities. We also explore the relevance of 
planetesimal size, comparing runs for a swarm of planetesimals 
having a radius of 10 and 100 km. Section [4] is devoted to dis- 
cussing our results in connection with other related studies and 
to summarising the main conclusions of our study. 

2. Outline of the overall model 

The description of the numerical model developed for the calcu- 
lation of giant planet formation and evolution was introduced in 
Benvenuto & Brunini (2005). Since then, several changes have 



Table 1. Characteristics of the Minimum Mass Solar Nebula at 
5.2 AU. 



Temperature 


133 K 


Solids surface density 


3.37 g crrr 2 


Gas volume density 


1.5 x 10" n gcm" 3 



been made to the code. In this section we will summarise the 
main improvements. 

2. 1 . The protoplanetary disc 

We characterise the protoplanetary disc with only three param- 
eters: the temperature profile, the solids surface density and the 
gas volume density, under the hypothesis of the model proposed 
by Hayashi ( fTMTT ). 

We consider a power law for the temperature profile 

r( fl ) = ar eff (i AU)r e >- 9 , (l) 

where a is the distance to the central star, T* is the effective 

' eft 

temperature of the central star, relative to that of the Sun (it is a 
dimensionless magnitude and in the cases considered in this ar- 
ticle T* s — 1), r e ff(l AU) is the effective temperature of the disc 
at 1 AU which was fixed at 280 K (the current temperature in the 
Solar System at that location), q is the temperature index (here 
q = 1 /2) and a is a factor that compensates for the fact that the 
central star was more luminous in the early epoch {a was elected 
to be 1.08). The snow line, a sn0 w, is located at T = 170 K, which 
corresponds in this model to a — 3.16 AU. 

The definition of the Minimum Mass Solar Nebula (MMSN) 
assumes a solids surface mass density profile that is a power law 

E(a) = 77 ice <r z (l AU) a^, (2) 

where we fixed <x z (1 AU) = 10 g ctrT 2 , p = 3/2 and 77; ce is the 
enhancement factor beyond the snow line, 

_ f 1 if a < a sn ow 
" ice I 4 if a > a snow . 

Although in this equation the assumption of a uniformly 
spread mass of solids is implicit, for the calculation of giant 
planet formation it is also necessary to adopt a planetesimal mass 
and radius (see the following sections). In this study we consider 
that the planetesimal bulk density, p m , is constant and equal to 
1 .5 g cm" 3 and that all planetesimals in the disc are of the same 
size. 

The gas volume density of the solar nebula is given by 

a-" 

p(a)=Og(l AU) — , (3) 

where <x g ( 1 AU) is the gas surface density at 1 AU taken to be 
2 x 10 3 g cm" 2 , b — 3/2 and H is the gas disc scale height, 

H(d) = 0.05 a 5/4 x (1.5 xl0 13 cm). (4) 

Note that a is always in AU but H must be in cm for p (a) to be 
in g cm" 3 giving the factor 1.5 x 10 13 cm. 

The main features of the MMSN at 5.2 AU are reported in 
Table □ 
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2.2. The core 

2.2.1 . The oligarchic growth of protoplanets 

Assuming that kilometre sized planetesimals can be formed in a 
protoplanetary disc, it is generally accepted that the first seeds 
for rocky protoplanets emerge through planetesimal accretion. 
The early stage in protoplanetary growth is the so-called "run- 
away growth" (Greenberg et al. 1978; Kokubo & Ida [1996). 
Runaway growth can be summarised as follows: random veloc- 
ities of large planetesimals are smaller than those of small plan- 
etesimals due to dynamical friction. This fact favours the gravita- 
tional focusing of large planetesimals which leads them to grow 
in a runaway fashion. These objects rapidly become more mas- 
sive than the rest of the planetesimals in the swarm and the first 
planetary embryos appear in the disc. During the runaway stage, 
both the growth rate of a protoplanet and the mass ratio of a pro- 
toplanet and planetesimals increase with time. 

Ida & Makino ( 119931 ) investigated the post-runaway evo- 
lution of planetesimals' eccentricities and inclinations due to 
scattering by a protoplanet. When runaway embryos become 
massive enough to affect planetesimals' random velocities, the 
growth regime switches to a slower one. Ida & Makino ( 1993) 
called this stage "the protoplanet-dominated stage", as proto- 
planets are now responsible for the larger random velocities of 
planetesimals which, in turn, decreases the growth rate of proto- 
planets. In their study of this regime, the authors performed 3D 
N-body simulations to investigate the evolution of eccentricities 
and inclinations of a system of equal-mass planetesimals and 
one protoplanet. Protoplanet-planetesimal and planetesimal- 
planetesimal interactions were both taken into account, but no 
accretion was considered. Ida & Makino ( 119931 ) found that plan- 
etesimals are strongly scattered by the protoplanet and part of 
the energy they acquire during the perturbation is distributed af- 
terwards among other planetesimals. The excited planetesimals 
determine a "heated region" around the protoplanet twice the 
width of the feeding zone of the protoplanet. In this region, the 
eccentricities e m and inclinations i m of planetesimals are highly 
perturbed and, as random velocities of planetesimals are well 
approximated by 



Vrel 



2\l/2 



)' /Z + <0 1/2 v k 



(5) 



the increase in e,„ and i m directly implies an increase in random 
velocities {\\ is the Keplerian velocity, {e 2 m ) 1 ^ 2 ((ij,) ) is the 
planetesimals RMS eccentricity (inclination)). 

Ida & Makino (119931 ) then proposed a two-step growth 
for the protoplanets. In the first stage, protoplanets grow in a 
runaway fashion in the usual sense: the stirring among plan- 
etesimals is dominated by planetesimals themselves and rela- 
tive velocities are low, favouring the rapid growth of the em- 
bryos. When a protoplanet is massive enough to perturb the 
surrounding planetesimals, the system enters the protoplanet- 
dominated stage. The numerical simulations show that this stage 
occurs in the dispersion-dominated regime and the relationship 
between mean planetesimals' eccentricities and inclinations is 
^m)'^ 2 /^™) 1 ^ 2 — 2- In this second stage, the growth rate of the 
protoplanet decreases as a consequence of the greater relative 
velocities between the protoplanet and planetesimals. However, 
the mass ratio of the protoplanet and planetesimals still increases 
with time. 

Ida & Makino (1993) also derived, using a semi-analytical 
model, the condition for the dominance of the protoplanet- 



planetesimal scattering over the planetesimal-planetesimal scat- 
tering: 



2 I M M > 2 m, 



(6) 



where M is the mass of the solid embryo, m is the effective plan- 
etesimal mass, S is the surface mass density of the planetesimal 
disc and 2m is defined as 



M 



-M 



2jraAa ' 



(7) 



with 2naAa the area of the "heated region". For a standard solar 
nebula, this condition can be translated into a relationship be- 
tween M and m which depends on the density profile of the disc, 
the semi-major axis a and the planetesimal mass. For the low- 
mass nebula model, the transition occurs when M/m ^50- 100. 
This condition was also corroborated by the authors through N- 
body simulations. Therefore, the second stage of the protoplanet 
growth begins when the mass of the protoplanet is still much 
smaller than the total mass of planetesimals in the feeding zone. 

Ida & Makino ( 119931 ) estimated the characteristic growth 
time (the mass-doubling time) for the protoplanet (assuming a 
constant solids surface density) : 



1 dM 
M dt 



oc M- 1/3 (e 2 m ). 



(8) 



To estimate (e 2 n ), they assumed an equilibrium state where the 
enhancement due to viscous stirring is compensated by the 
damping due to the nebular gas drag. As a result, in the first stage 
where the dominant perturbers are the planetesimals, {e 2 /) oc 
m 8/i5 ( note tn at j t does not depend on the mass of the proto- 
planet), while in the second stage {ef n ) oc m ' M 2 ' 3 . Hence, in 
the first stage T grow oc M -1 ^ 3 so protoplanets grow very fast and 
the ratio M/m =s 50 - 100 is reached in a negligible time. In 
the second stage, the dominant perturber is the protoplanet and 
Tgrow 00 M 1 ^, thus the growth of the protoplanet gradually be- 
comes slower as its mass increases. This type of time depen- 
dence is characteristic of the orderly growth (the presumed final 
growth regime for the terrestrial planets). However, the growth 
regime in the protoplanet-dominated stage is much shorter than 
the orderly growth because protoplanets grow by accretion of 
planetesimals and not through collisions of comparable sized 
bodies where dynamical friction is absent. The three growth 
regimes (the runaway regime, the 2-step regime and the orderly 
regime) are schematically illustrated in Fig. 1 1 of Ida & Makino 
(TT9931 . 

Kokubo & Ida (1998) performed 3D N-body calculations 
to study the post-runaway accretion scenario. Their simula- 
tions start with two protoplanets that grow by planetesimal ac- 
cretion. In the first stage of accretion, protoplanets grow in a 
runaway fashion. Protoplanetary runaway growth breaks down 
when M ~ 50 m, as predicted by the semi-analytical theory of 
Ida & Makino ( I1993I ). Protoplanets subsequently grow keeping 
a typical orbital separation of 107?h, where Ru is the Hill radius 
of a protoplanet, 



R« = 



M 
3MT 



1/3 



(9) 



with a the distance to the central star. The orbital repulsion is a 
consequence of the coupling effect of scattering between large 
bodies and dynamical friction from the swarm of planetesimals. 
Protoplanets continue their growth keeping their mass ratio close 
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to unity, but the mass ratio between protoplanets and planetes- 
imals continues to increase with time. Thus, protoplanets grow 
oligarchically, as no substantial accretion between planetesimals 
themselves is found. The mass distribution then becomes bi- 
modal: a small number of protoplanetary embryos and a large 
number of small planetesimals shape the protoplanetary disc. 
Kokubo & Ida d 19981 1 coined this growth regime "oligarchic 
growth", 'in the sense that not only one but several protoplan- 
ets dominate the planetesimal system' . 

Kokubo & Ida (2000 ) investigated, through 3D N-body sim- 
ulations, the growth from planetesimals to protoplanets includ- 
ing the effect of the nebular gas drag. They confirmed the exis- 
tence of an initial runaway phase in protoplanetary growth and 
a second stage of oligarchic growth of protoplanets, in the same 
sense as in Kokubo & Ida ( 119981 1. Also, when analysing the evo- 
lution of the RMS eccentricity and inclination of the planetes- 
imal system during the oligarchic growth stage they found that 
their results agree with those predicted by the semi-analytical 
theory of Ida & Makino ( 1993 ) for their second stage. 



2.2.2. The growth of the core 

In the present version of the code we have introduced a time- 
dependent planetesimal accretion rate. In this kind of calculation 
most author s (Pollack et al l 19961 Hubickyj et al. 120051 Alibert 
et al. 120051 1 usually prescribe that obtained by Greenzweig & 
Lissauer (1992) which assumes a rapid growth regime for the 
core. Instead, we adopt that corresponding to the oligarchic 
growth of Ida & Makino ( 1993 ); a slower accretion rate that still 
has not been explored with a self-consistent code for giant planet 
formation. 

The condition for the dominance of the oligarchic growth 
over the (previous) runaway growth of a protoplanet, M/m ^ 
50-100, was derived semi-analytically by Ida & Makino ( 11993I ). 
In all the cases of interest for this study, the cross-over mass is 
very low — i.e. some orders of magnitude below the Earth mass 
(Thommes et al. 120031) . Due to the initial runaway regime, the 
cross-over mass is reached in a negligible time and thus, the 
formation time of a giant planet's core is almost entirely regu- 
lated by the oligarchic growth. For this reason, we prescribe the 
oligarchic growth for the core since the very beginning of our 
simulations. 

In the dispersion-dominated regime, a solid embryo growth 
rate is well described by the particle-in-a-box approximation 
(Safranov. [T969l ). 



dM c 
dt 



(10) 



where M c is the mass of the solid protoplanet (in our case, the 
core of the giant planet), h is the planetesimals' disc scale height, 
R e s is the effective capture radius and F is a factor introduced 
to compensate for the underestimation of the accretion rate by a 
two-body algorithm when considering the velocity dispersion of 
a population of planetesimals modelled by a single eccentricity 
and inclination equal to the RMS values. F is estimated to be 
~ 3 (Greenzweig & Lissauer [1992]). 

Due to gravitational focusing, the effective capture radius of 
a protoplanet, R e ff, is larger than its geometrical radius 



d2 r>2 



1 + 



2\ 



(ID 



with R c the geometrical radius of the solid embryo, v esc the es- 
cape velocity from its surface and v re i the relative velocity be- 
tween the protoplanet and planetesimals, 



V'rel 



Ve 2 + f a Q k , 



(12) 



where, hereafter, e = (e^) 112 (i = (if n } 1 ^ 2 ) is planetesimals RMS 
eccentricity (inclination) with respect to the disc (<?, i « 1) and 
is the Keplerian angular velocity. We apply the approxima- 
tions i ^ e/2 and h ^ ai. Following Thommes et al. (2003), we 
adopt for e the equilibrium expression that is deduced for the 
case when gravitational perturbations due to the protoplanet are 
balanced by dissipation due to gas drag, 



(13) 



where M is the protoplanet mass (here we assume M to be the 
total mass of the proto-giant planet, meaning the core mass plus 
the envelope mass), p m is the planetesimal bulk density, p is the 
gas volume density of the protoplanetary disc, Cd is the drag co- 
efficient (dimensionless and of the order of 1), M* is the mass of 
the central star and /3 is the width of the "heated region" in units 
of the Hill radius (considering that potentially other embryos are 
growing as well, /3 is of the order of 10). For calculations that do 
not assume equilibrium values for the eccentricity and the incli- 
nation, we refer the reader to Chambers (2006). His results show 
that using equilibrium values for e and i is an acceptable approx- 
imation when considering embryos of m ^ 10~ 2 M ffi , consistent 
with the initial core mass of all our simulations (see Sect. [3). 

However, when the solid embryo is massive enough to grav- 
itationally bind gas from the surrounding nebula, the presence 
of this envelope should be considered when calculating the ef- 
fective capture radius, R e ff. The gaseous envelope modifies the 
trajectory of incoming planetesimals, as they are affected by the 
gas drag that enlarges the capture radius of the protoplanet. In 
addition to the gravitational focusing, a "viscous focusing" due 
to gas drag should be considered. 

The effective radius R e ff, in the form stated in Eq. ( fTTT i dom- 
inates when the mass of bound gas is negligible. But when the 
embryo acquires enough gas to form a thin atmosphere, its effec- 
tive radius becomes larger and it separates from that calculated 
previously. 

To calculate the effective radius of the protoplanet in the 
presence of gas around the solid embryo we take into account 
the action of gravity and gas drag on the incoming planetesimals. 
Consider one planetesimal entering the protoplanet atmosphere 
with a velocity v re i (Eq. (fT2l). Its equation of motion results from 
the action of gravity 



fc 



GM r m, 



(14) 



together with the action of gas drag. In Eq. (|14t . r is the radial 
coordinate from the centre of the protoplanet and M r is the mass 
contained within r. For the gas drag force acting on a spherical 
body of radius r m travelling with velocity v, we adopt the Stokes 
law 



/d = --C D nrlp g v 2 v, 



(15) 



with p g the envelope density and Cd = 1 (Adachi et al. 119761 ). 
Starting at the external radius of the protoplanet (see Sect. 12.41 
Eq. (f27j) for its definition), numerical integrations of the equation 
of motion are performed with an adaptive Runge Kutta 4 routine 
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Fig. 1. Comparison between the protoplanet effective radius for 
planetesimal capture. In solid line, R* ff is calculated using Eq. 
(fTTT i and in dashed-dotted line R e g is calculated including the 
gas drag of the atmosphere. These results were obtained with a 
full evolutionary numerical simulation. The protoplanetary disc 
surface density corresponds to a 6 MMSN. The protoplanet is 
located at 5.2 AU and the radius of incoming planetesimals is 
100 km. 



(Press et al. |1992| l. The initial impact parameter considered for 
the integration of the trajectory is the effective radius in the ab- 
sence of gas (Eq. ( fTTT ). hereafter 7?* ff ). Note that this is the lowest 
possible effective radius of the protoplanet. We calculated the 
trajectory of the planetesimal until it completes one close orbit 
inside the Hill sphere of the protoplanet (so it is always gravita- 
tionally dominated by the protoplanet). At this point, the plan- 
etesimal is considered captured and will definitely be accreted 
after a certain time. This procedure is repeated for increasing im- 
pact parameters (the relative difference between two consecutive 
trial impact parameters being 5 x 10~ 4 ) until the resulting orbit 
is no longer inside the Hill sphere. In other words, planetesimals 
are considered to be captured when the energy lost by gas drag 
allows them to complete one close orbit around the core and this 
is completely inside the protoplanet. The largest impact parame- 
ter for which this condition is fulfilled is adopted as the effective 
radius for that model. The solids accretion rate is still Eq. dlOt . 
using this new definition for R e g, The mass of accreted planetes- 
imals is added to the mass of the core. It is worth mentioning 
that according to Inaba & Ikoma (2003 ), a detailed calculation 
of planetesimal ablation does not have a significant influence on 
the estimation of /? e ff ■ 

Note that our criterion for the capture of planetesimals and 
the determination of R e ff is different to that employed by Pollack 
et al. ( 1 19961 1 where R e ff is taken as the periapsis altitude of the 
outermost bound orbit. 

We performed a full evolutionary calculation of the forma- 
tion of a giant planet in circular orbit around the Sun, with 
a = 5.2 AU, and immersed in a disc for which solids and gas 
densities are increased by a factor of 6 with respect to the MMSN 
(see Sect. |3.21 >. The effect of gas drag on the protoplanet effective 
radius is very important, as seen in Fig. Q] The effective radius 
calculated including the gas drag separates from the effective ra- 
dius calculated according to Eq. ( fTTT i. R* w , when the elapsed time 
is t ~ 4 My and the radius of the core is R c ~ IQ 9 cm. This cor- 
responds to an envelope mass of ~ 5 x 10~ 4 M ffi and a core mass 
of ~ 2 M e . Note that by the end of the formation process, R e g is 



more than 10 times larger than R* s . Moreover, it enters quadrati- 
cally in the core accretion rate (Eq. (fTOl)), so it will impact on the 
protoplanet formation timescale and on the final core mass (see 
Sect. gj. 

The core growth rate (Eq. ( fTOl )) depends also on the surface 
mass density of planetesimals, S, a quantity that is neither con- 
stant in time nor in space. E is a function of the distance to the 
central star (see Sect. 12.11 ) and also depends on the disc evolu- 
tion. In this paper we consider only one effect for the time de- 
pendence of S: the accretion by the protoplanet. There are other 
important effects to be taken into account for a more realistic 
calculation of the planetesimal accretion rate, for example, plan- 
etesimal migration in the protoplanetary disc due to the nebu- 
lar gas drag (Thommes et al. 2003 ) and planetary perturbations, 
for example, planetesimal ejection away from the feeding zone 
of the protoplanet (Alibert et al. 2005 ). However, this paper is 
mainly focused on the coupling of an oligarchic growth regime 
for the core with a full evolutionary calculation of the formation 
of a giant planet, so in this first exploration we shall neglect other 
effects that may contribute to the time evolution of S. 

When considering the time variation of 2 we assume that the 
growing protoplanet is able to capture planetesimals only from 
its feeding zone, defined as an annulus of 4 Hill radii (Eq. ((9])) 
at either side of its orbit. The surface density inside the feeding 
zone at a certain time t is considered to be uniform, changing 
only due to the accretion by the protoplanet 



Y = <£> (?) = <£> 



M c (f) 



(16) 



where a mt (t), a out (t) are the inner and outer boundaries of the 
feeding zone and (S) is the average of the initial surface den- 
sity in that region. Adopting this expression for E tacitly implies 
that the accretion process is slow enough to guarantee a uniform 
distribution of planetesimals inside the feeding zone. 

2.3. Interaction between incoming planetesimals and the 
protoplanetary envelope 

On their trajectory to the core, accreted planetesimals interact 
with the gaseous envelope of the protoplanet exchanging energy. 
In this study, we implemented a simple model to take into ac- 
count this effect (for more detailed models see Podolak et al. 
[19881 Pollack et al. [T9961 Alibert et al. [20031 . To simplify the 
situation, we assume planetesimals approaching the protoplanet 
as coming from infinity and entering the envelope with velocity 
v le i (Eq. (fT2li). describing afterwards a straight trajectory to the 
core. This may seem to contradict the fact that we calculate the 
orbits of planetesimals to define R e s, but as mentioned earlier we 
stop orbit calculation when planetesimals complete the first rev- 
olution inside the Hill sphere. In the future we will incorporate 
a complete calculation of the trajectories to the core, together 
with a self-consistent deposition of planetesimal energy in the 
gaseous layers of the envelope. 

In our simple model of a straight trajectory to the core, the 
energy of one incoming planetesimal of mass m at the top of the 
atmosphere is 



1 9 

E = —m v„ 



(17) 



Gas drag and gravity acceleration are considered here as the two 
main responsible forces for changing the energy of planetesimals 
throughout the envelope. The drag force acting on a spherical 
planetesimal of radius r m and moving with velocity v is given 
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by Eq. (JT3J. The gravitational force is calculated as usual with 
Eq. CEl. 

Both forces act simultaneously and modify the planetesimal 
velocity. The total force acting on a planetesimal is 



/ = /d+/c 



dv 
-m — r 
dt 



and 

dv dv 
dt dr 7 

then 

dv Co7irf n p g v GM r 



dr 



2m 



(18) 



(19) 



(20) 



The kinetic energy lost by a planetesimal due to gas drag 
£k,m is transformed into heat, gained by the envelope's shells, 



dE 
dr 



dE\ 



k.m 



1 



dr 



C D nr 2 m p $ v 2 . 



(21) 



These equations correspond to the energy exchange of one plan- 
etesimal with the surrounding gas. When considering all incom- 
ing planetesimals per unit of time, M c , and transforming the dif- 
ferential equations into difference equations, the total variation 
of energy AE, of the envelope's shell i in the time interval At is 



1 9 M c , 

AEi = -C D nr z mPg —vf An At. 

2 m 



(22) 



where v, (the velocity of planetesimals when entering an enve- 
lope's shell z) is obtained from a discretisation of Eq. ( 1201 ). 



Vi = V;-i + 



C D xrlp g v,-_i GM n _ ]/2 ' 



r i-l/2 V '-i) 



2m 



An-i. 



(23) 



Finally, when planetesimals reach the core, all their remain- 
ing kinetic energy is deposited in the adjacent gas shell. 

2.4. The gaseous envelope 

The general procedure employed in this study to numerically 
solve the formation of a giant planet was introduced in a pre- 
vious paper (Benvenuto & Brunini 2005). In the above sections 
we mentioned some of the main improvements introduced to the 
model presented in that work. For the sake of completeness, we 
summarise here the fundamental ideas concerning the solution of 
the envelope's structure. The full differential equations of giant 
planet formation are solved with an adaptation of a Henyey-type 
code. Radiative and convective transport are considered, em- 
ploying the Schwarzshild criterion for the onset of convection. 
The adiabatic gradient is adopted for the temperature gradient in 
the latter case. 

The boundary conditions remain the same. As inner bound- 
ary conditions, we consider the core density to be constant, 
p c = 3 g cm 4 . The mass of the core at time t is calculated as 



M c = -np c Ri(t). 

The luminosity on the surface of the core is 
L r (M c ) = 0, 



(24) 



(25) 




Time [10 years] 

Fig. 2. Comparison case. We plot the mass evolution of the core 
(solid line), of the envelope (dash-dotted line) and of the total 
mass (dashed line) of the protoplanet. The input parameters of 
this run are the same as those of case J3 of Pollack et al. ( 1996). 
This figure should be compared to Fig. 2b of the cited article. 
See the text for details. 



and the velocity is 

Mc 



v(M c ) 



(26) 



For the outer boundary conditions, the planetary radius is 
defined as 



R = min[/f acc ,/? H ], 

where the accretion radius, /? a cc, is given by 



Rj 



GM 



(27) 



(28) 



c being the sound velocity. Ru is the Hill radius (Eq. ©). The 
external boundary conditions for the temperature and density of 
the envelope are those corresponding to the nebular gas, T and 
p (see Sect. 12. 11 1. For further details on this aspect, we refer the 
reader to Benvenuto & Brunini d20051 l. 

Regarding the equation of state (EOS), we employ results 
presented in Saumon et al. (1995). For the grain opacities, we 
use the tables from Pollack et al. ( 1985). For temperatures above 
10 3 K we consider Alexander & Ferguson d 19941 ) molecular 
opacities which are available up to T < 10 4 K and for higher 
temperatures we consider opacities by Rogers & Iglesias d!992| l. 

3. Results 

3. 1 . Comparison with previous studies 

To our knowledge, no self-consistent calculation of giant planet 
formation with a core growing according to the oligarchic model 
has been performed to date. To show the impact of the core ac- 
cretion rate on the whole formation of a planet we performed 
a simulation to be compared with previous results obtained by 
authors adopting another core accretion model. 

For comparison we selected one of the simulations per- 
formed by Pollack et al. ( |1996l l, their case J3 for the formation 
of a Jupiter-like object. These authors adopted for the growth of 
the core the model of Lissauer (1987 ), which assumes a more 
rapid growth regime than the model of Ida & Makino (1993) 
adopted here. For this run, the protoplanetary disc at 5.2 AU and 
the accreted planetesimals are characterised as follows: 
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Fig. 3. For the same comparison case as in Fig. [2] solid (dashed- 
dotted) line shows the time evolution of the logarithm of the core 
(envelope) growth rate. 



- initial solids surface density 2 = 15 g ctrT 2 , 

- nebula volume density p = 5 x 10~ n g ctrT 3 

- nebula temperature T — 150 K, 

- planetesimal bulk density p,„ = 1.39 g cm" 3 , 

- planetesimal radius r m = 100 km. 




Fig. 4. The impact of the selected model for the growth of the 
core. The solid line shows the ratio of eccentricities e/e*, where e 
is calculated according to Eq. ( fT3l and <?„ according to the model 
adopted by Pollack et al. dl996l l. The dashed-dotted line plots 
the ratio of inclinations ///*. While in our model i depends on 
e (i - e/2), in Pollack et al.(1996) i, remains constant during 
the entire simulation. For the X axis we selected the Hill radius 
in order to get rid of the different formation timescales of both 
simulations. 



The initial mass of the core, M c , is 0.1 M ffi and the core density, 
p c , is 3.2 g cm 4 . 

The resulting core and envelope mass evolution of our run 
is depicted in Fig.|2j which should be compared with Fig. 2b of 
Pollack et al. ( 1996). We use the same mass range in the Y axis 
to favour the comparison, although the run was completed suc- 
cessfully to the end, that is, until the total mass of the planet, M, 
is that of Jupiter (318 M ffl ). From Fig . [2] it can be seen that in our 
calculation it takes 19 My for the mass of the envelope to equal 
the mass of the core, while in Pollack et al. i 19961 1 the cross- 
over point (M c = M env ) is reached in only 1.51 My. As is imme- 
diately obvious, the time difference between both calculations is 
of approximately one order of magnitude. This is a significant 
difference if we keep in mind that the estimated lifetime of the 
solar nebula is less than 10 My. However, it is worth pointing out 
that the cross-over mass in both simulations is almost the same 
(~ 29 M ffi ). For the sake of completeness, we mention that the 
whole formation time is 20.6 My and the final mass of the core is 
42 M ffi . From the comparison of both figures, it can also be seen 
that we do not find three phases in the formation of the planet as 
Pollack et al. ( 1996 ) found in their simulations. The slow growth 
of the core guarantees a smooth variation of the slope of the M c 
curve. No different growth regimes in the mass of the core or 
in the mass of the envelope can be distinguished before the run- 
away growth of the envelope. Moreover, in Fig. [3] we show in a 
logarithmic scale, the solids and the gas accretion rate. Pollack 
et al. ( 119961 ) obtained a different behaviour of the accretion rates 
which allowed them to define a short phase 1 in the growth of the 
protoplanet which ends when dM c /dt = dM env /dt, and a longer 
phase 2 which ends when M c = M env . Phase 2 is characterised 
by a constant proportionality between dM c /dt and dM env /dt. 
Phase 3 corresponds to the usual runaway growth of the enve- 
lope. From Fig. [3] we can see that a distinction between phase 1 
and a phase 2 cannot be inferred from our simulation. The time 
derivative of the mass of the envelope is a monotonously increas- 
ing function, with no flat slope after it becomes larger than the 
core growth rate. 



Our physical model is in essence the same as that of Pollack 
et al. (119961 ), except for the selection of the core accretion model 
and the treatment of the interaction between incoming planetesi- 
mals and the gaseous envelope of the protoplanet. Regarding the 
latter, as explained in Sect. 12.31 we consider a simple model for 
the energy deposition, while Pollack et al. (1996) developed a 
more complex model where the complete trajectories to the core 
of accreted planetesimals are calculated and planetesimal abla- 
tion is also taken into account. However, this can not be the main 
cause of such an important difference between the two formation 
timescales. Although this effect may be important in determin- 
ing the internal structure of the protoplanet and may also modify 
the formation timescale, for this to happen a significant amount 
of gas bound to the core is needed. In the present simulation, 
it takes 10 My to form a core of 10 M ffl , and the mass of the 
envelope at that stage is only 0.1 M Q . Hence, the main cause 
of the significant discrepancies between the timescales, must be 
the growth model for the core. We adopt an oligarchic growth for 
the core (with an enhanced effective radius due to the gas drag 
of the envelope), while Pollack et al. ( 1996 ) prescribed the more 
rapid accretion model of Lissauer ( 119871 ), also modified by the 
enhancement caused by the gas drag. According to their Eq. (1), 
the core accretion rate can be written as: 

dM , 

— = nR 1 I.D.F f , (29) 
dt 8 

where R is the effective radius, D. is the Keplerian angular veloc- 
ity and F g is the gravitational enhancement factor. F g depends 
on the RMS values of the reduced eccentricities and inclinations, 
and on the reduced effective radius, 

a . a . _ R 

K H K H Kr 

They calculate the numerical value s of F a using the formulae ob- 
tained by Greenzweig & Lissauer ( 1992). Their model assumes 
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that planetesimals' inclinations depend only on planetesimal- 
planetesimal interactions and adopt for i^ m the formula: 



where v e , m is the escape velocity from the surface of a planetes- 
imal. Clearly, the inclination z'„ = i^,„ Rn/a remains constant 
during the entire calculation while decreases with time. On 
the other hand, they assume that eccentricities are controlled by 
both planetesimals and protoplanet stirring, and prescribe for the 
reduced eccentricity the equation 

e h ,m = max(2z' h , m ,2). (31) 

This means that if <?h, m = 2z'h )BJ the protoplanet is growing ac- 
cording to the runaway regime, as e and i would be indepen- 
dent of the mass of the protoplanet, while if eh,m = 2 the ec- 
centricity of planetesimals would be affected by the presence of 
the protoplanet, this condition corresponding to the protoplanet- 
planetesimal scattering in the shear-dominated regime. The 
numerical simulations of Ida & Makino (1993) showed that 
the protoplanet-planetesimal scattering in the shear-dominated 
regime lasts for only a few thousand years, after which plan- 
etesimals are strongly perturbed by the protoplanet, and the 
post-runaway regime can then be considered as a dispersion- 
dominated regime. In the dispersion-dominated regime eccen- 
tricities and inclinations satisfy that eh,m/*iyn — 2 and e\ m > 2. 
Hence, in a model where the protoplanet-planetesimal scattering 
is in the shear-dominated regime, eccentricities and inclinations 
of planetesimals in the vicinity of the protoplanet remain low. 
This leads to an accretion scenario which is much faster than 
that corresponding to the oligarchic regime. 

To show the difference between eccentricities and inclina- 
tions in both models throughout the formation of the planet, we 
included in our simulation of the case J3, the calculation of i^,„ 
and eh,m according to Eqs. (f30b and OTb . The ratio of the eccen- 
tricity and inclination between the two models is shown in Fig. 
|U We chose the X axis to be the Hill radius of the protoplanet, as 
it is completely independent of time and the comparison is then 
straightforward. For this case, the initial value of is lower 
than 1 and, as it is a decreasing function of the mass of the pro- 
toplanet, eh,m always equals 2. Consequently, our eccentricities 
are larger than those in Pollack et al. ( 1996 ) by a factor of 4, as 
can be seen from Fig. [4] On the other hand, our inclinations are 
much higher (their inclination is constant throughout the whole 
formation of the planet, z'„ 0.0039), which impacts not only in 
the relative velocities but also in the scale height of the planetes- 
imal disc which is inversely proportional to the core accretion 
rate. As a consequence, due to our larger eccentricities and in- 
clinations, we obtain a much lower accretion rate and a longer 
formation timescale. 

Pollack et al. (1996 ) adopt for the definition of the bound- 
aries of the feeding zone that the radial distance on either side 

of the orbit is given by Aa = ^j\2 + e^ m Rn. According to our 

simulation, e^„, is always equal to 2, which leads to Aa = 4-Ru, 
the same value we adopted in our calculations. 

3.2. Results of full calculations of giant planet formation with 
an oligarchic growth for the core 

All the simulations presented in this section were started with 
an embryo of 0.005 M ffi , revolving in a circular orbit around the 
Sun. The semi-major axis of the orbit is that corresponding to 



Jupiter, 5.2 AU. The mass of the seed was chosen to guarantee 
that the embryo is undergoing an oligarchic growth and that it is 
also able to bind a tiny atmosphere, the latter in order to allow the 
code to converge. The density of the core will be held constant 
during the entire formation process, p cole = 3 g cm 4 . The dis- 
crepancy between core density and planetesimals bulk density 
(p m — 1 .5 g cm 4 ) is due to the progressive high pressure the 
core is subject to. We arbitrarily stopped calculations when the 
total mass of the protoplanet was that of Jupiter. This assumption 
tacitly implies that the nebula can always provide the necessary 
amount of gas for the formation of the planet. 

The first set of simulations is intended to show the conse- 
quences of the envelope's drag on the oligarchic growth of the 
core (Eq. (fTOli). To this end, we performed several runs for vari- 
ous protoplanetary disc surface densities, including and exclud- 
ing the effect of the gas drag in the calculation of R e s. The plan- 
etesimal radius is assumed to be constant, fixed for these exam- 
ples at 100 km (no size distribution or collisions between plan- 
etesimals were taken into account). As no planetesimal migra- 
tion, planetesimal ejection by the protoplanet, perturbations by 
other embryos, etc., were considered, the variation of the solid 
surface density of the feeding zone is only due to planet accre- 
tion. Since no mechanisms for the dissipation of the gaseous 
component of the nebula were modelled, the variation of the vol- 
ume gas density of the feeding zone is due only to the formation 
of the envelope. We note that this is an idealised situation where 
the protoplanet probably has the largest amount of available ma- 
terial to feed itself. 

We found that, in spite of the favourable conditions for mass 
accretion of our model, for protoplanetary discs with densities 
lower than that of a 6 MMSN, the formation process could not 
be completed according to the timescales imposed by the ob- 
servations of circumstellar discs (lower than 10 7 yr) in either of 
the cases considered for the calculation of R e $. The results for a 
disc of 6 MMSN are depicted in Fig. [5] The upper panel (Fig. [5] 
a) shows the evolution of the core and envelope mass when ig- 
noring for the protoplanet capture effective radius the presence 
of the atmosphere (R* eff calculated according to Eq. (fTTTi). The 
complete formation of a Jupiter-mass object takes a bit over 17 
My and the final mass of the core is 28 M ffi (note that, in our 
model, all accreted solids are deposited onto the core). However, 
when including the atmospheric gas drag, the timescale turns out 
to be 12 My (still over the limiting 10 My) while the mass of the 
core increases to 60 M e (see Fig.|5]b). This means that, in this 
case, the effect of the gas drag of the atmosphere reduces the for- 
mation time in about 30%, but also affects the final mass of the 
core, which is increased a factor of 2. 

Although a core of ~ 10 M e (currently an acceptable value 
for Jupiter's core mass, Saumon & Guillot 2004) is formed in 
~ 6.5 My in the second simulation, the runaway collapse of the 
gaseous envelope occurs ~ 5 My later. As seen from Eq. ( TTOb . 
the solids accretion rate is directly proportional to the solids sur- 
face density and the fact that the feeding zone is far from being 
depleted when this mass is achieved (see Fig. [6j, allows the pro- 
toplanet continue the accretion of solids. In fact, as seen from 
Fig. [7] the peak of the core accretion rate takes place for t ~ 8 
My, when the mass of the core is approximately 20 M®. This 
large number of incoming planetesimals contribute to the gas 
pressure supporting the envelope via their gravitational energy, 
and thus prevent faster accretion of gas from the nebula and as 
such permit the formation of a massive core. 

The same behaviour is found when increasing the surface 
density of the disc. In Table [2] we show, in round numbers, the 
formation times and the final masses of protoplanetary cores for 
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Fig. 5. For a 6 MMSN protoplanetary disc, panel a shows the 
core (solid line) and envelope (dashed-dotted line) mass evolu- 
tion of a Jupiter-like forming protoplanet at 5.2 AU, considering 
that the core grows according to the oligarchic growth regime, 
with R*~ calculated with Eq. (fTTT i. Panel b shows the same sim- 
ulation but including the gas drag effect of the protoplanetary 
atmosphere in the calculation of R e g. Note that for a clear visu- 
alisation the Y axis changes from a linear scale to a logarithmic 
scale. Panel a (b) is in a linear scale up to 40 M e (100 M e ), 
afterwards the scale is logarithmic. For both runs, planetesimal 
radius was set equal to 100 km. 
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Fig. 7. This figure shows the evolution of the core accretion rate 
for the simulation that includes atmospheric gas drag. Note that 
the accretion rate is still an increasing function when the mass 
of the core reaches 10 M ffi (corresponding to ~ 6.5 My; see also 
Fig- 0b), preventing the rapid accretion of the envelope mass. 

Table 2. Comparison of formation times and core masses for 
several disc densities. Planetesimal radius is 100 km. 



Disc density 


Without gas drag 


With gas drag 


[MMSN] 


Formation 


Core 


Formation 
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time [My] 
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Table 3. Comparison of formation times and core masses for 
several disc densities with a planetesimal radius of 10 km. 



Disc density 
[MMSN] 



Without gas drag 



With gas drag 



Formation Core Formation Core 
time [My] mass time [My] mass 
[M e ] [M e ] 
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Fig. 6. This panel depicts the evolution of the solids surface den- 
sity for the simulation that includes the atmospheric gas drag 
(corresponding to Fig.[5]b). A core of ~ 10 M e is achieved for 
t ~ 6.5 My, when 2 is still very high. 
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discs that range from 6 MMSN to 10 MMSN and planetesimals 
having a radius of 100 km. As expected, the higher the density 
of the disc, the shorter the formation time but the larger the final 
core's mass (see also Pollack et al. II 9961) . The timescale reduc- 
tion, when the gas drag from the atmosphere is included, ranges 
from ~ 30% for 6 MMSN to ~ 40% for 10 MMSN. The maxi- 
mum timescale reduction for these cases is then less than a factor 
of 2. The mass of the core, on the other hand, monotonously in- 
creases fractionally, from a factor of ~ 2.15 for 6 MMSN to a 
factor of ~ 2.35 for 10 MMSN. Qualitatively, all curves look 
similar to that shown in Fig. [5] and it is not worth showing them 
here. 

When the same comparative simulations are made for plan- 
etesimals having a radius of 10 km (this means, simulations that 
consider both cases for the effective radius of the protoplanet, 
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Fig. 8. For a protoplanetary disc of 6 MMSN and planetesimals 
having a radius of 10 km, upper panel (a) shows the evolution of 
the protoplanet mass (the core mass with the solid line, the en- 
velope mass with the dashed-dotted line) when no atmospheric 
gas drag is considered. The lower panel (b) illustrates the same 
process but when is calculated including the gas drag force 
acting on incoming planetesimals. Again, the Y axis is split in 
two different scales: first we adopt a linear scale and then a log- 
arithmic one. The same range as in Fig. (0 is adopted. 



R* efi and /? e ffX the timescale reduction ranges from ~ 38% for 6 
MMSN to ~ 56% for 10 MMSN, and again we find that the mass 
of the core is doubled (see Table [3j. For planetesimals of this 
size, and according to the results presented in Table[3] one could 
speculate about reducing the density of the disc (to less than 
that corresponding to 6 MMSN) and still finding a formation 
timescale below 10 7 yr (when considering the enhanced capture 
radius by gas drag). However, to present our results clearly, we 
find it more relevant to compare the several sets of simulations 
for the same nebula conditions. 

For the sake of completeness, Fig. [8] shows the evolution of 
the core mass and the envelope mass for a 6 MMSN disc where 
the planetesimal radius is 10 km. When gas drag from the en- 
velope of the protoplanet is ignored, the formation time is ~ 9 
My and the core mass is ~ 34 M ffl . However, when this effect is 
included in the calculations, the total formation time is reduced 
to ~ 6 My, with a final super-massive core of 69 M®. 

Although our core masses are in all cases completely out 
of range for present models of Jupiter total mass of solids, Mz 
(recent estimations are M c ^ 1 1 M e and 1 M e ^ M Q + M z ^ 
39 M e , Saumon & Guillot 2004), the main objective of this work 
is not directed to fit our results to these values, but on stressing 



the results obtained when adopting the oligarchic growth for the 
core. From our calculations it is clear that including the enhance- 
ment of the accretion cross-section due to gas drag, helps to 
reduce the formation timescale when considering an oligarchic 
model, although surface densities higher than that of the MMSN 
are still needed. The main drawback of these massive discs is 
that they lead to the formation of super-massive cores. However, 
we have to keep in mind that no mechanisms of planetesimal re- 
moval from the feeding zone of the protoplanet are taken into 
account, except for the depletion due to planetesimal accretion 
onto the protoplanet. 

Thommes et al. ( 120031 ) calculated analytically and numeri- 
cally the evolution of solid protoplanetary masses growing oli- 
garchical for a wide range of semi-major axes and for a pro- 
toplanetary disc of 1 MMSN and 10 MMSN (their solids sur- 
face density for the MMSN is slightly lower than ours). For a 10 
MMSN and considering planetesimals having a radius of 10 km, 
they obtain a protoplanet of 60 M©, a result we were also able 
to reproduce in the absence of gas accretion. However, when 
they include in their simulations the planetesimal migration in 
the protoplanetary disc due to nebular gas drag, the resulting 
protoplanet is reduced to ~ 10 M e . Planetesimal migration due 
to gas drag is a relevant ingredient which should be considered 
in future calculations for an accurate estimation of the mass of 
the core and of the formation timescale. Another important ef- 
fect that should be taken into account in future calculations is 
the ejection of planetesimals out of the feeding zone. This could 
probably also help in solving the problem of having a Jupiter- 
like planet with such large core masses. Since ejection occurs 
mainly for large mass planets, this would not necessarily slow 
down the accretion of planetesimals at the beginning, and may 
not lead to a too long formation time of the planet. 

The previous paragraphs focus on the consequences of the 
enhancement of the effective capture cross-section due to the 
envelope's gas drag. We now analyse the importance of plan- 
etesimal size. For a clear visualisation we summarise in Table 
[4] the earlier calculations with the enhanced effective radius for 
planetesimals having radii of 10 and 100 km. Planetesimal mass 
appears in Eq. ( fT3] l, which determines the relative velocity of 
planetesimals with respect to the protoplanet through Eq. ( fT2l . 
Smaller planetesimals experience stronger damping of random 
velocities due to nebular gas and form a thinner disc (i st j e 
and h =s a i) which in turn favours accretion (Thommes et al. 
120031 1. Moreover, smaller planetesimals suffer from a stronger 
deceleration due to atmospheric gas drag, increasing the impact 
parameter. Both effects combined favour the formation process, 
reducing the formation timescale. 

If we analyse and compare results without including the gas 
drag of the envelope for a fixed nebula density (first and second 
Cols, of Tables [2] and [3j, we find that for planetesimal radius 
of 10 km the whole formation process is reduced in an approxi- 
mately constant factor of 2 when compared to planetesimals hav- 
ing a radius of 100 km, while the core mass is increased by 25%. 
This value is independent of the protoplanetary disc conditions, 
at least in the cases considered here. When the gas drag is in- 
cluded in the calculations this behaviour changes (see Table 0J. 
Results show that the formation time reduces from a factor of ~ 2 
for a protoplanetary disc of 6 MMSN to a factor of ~ 3 in the 
case of 10 MMSN, when comparing the runs for the same pro- 
toplanetary disc conditions but for different planetesimal size. 
The factor of the timescale reduction depends then, on the sur- 
face density of the protoplanetary disc. The mass of the core is 
less affected than in the previous case, resulting for planetesi- 
mals having a radius of 10 km, as being on average, 15% more 
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Table 4. Comparison of formation times and core masses for 
several disc densities, and for two different planetesimal sizes: 
100 and 10 km. All cases were calculated including atmospheric 
gas drag. 



Disc Planetesimal Planetesimal 
density radius: 100 km radius: 10 km 



[MMSN] 


Formation 


Core 


Formation 


Core 




time 


mass 


time 


mass 




[My] 


[M e ] 


[My] 


[M ffi ] 


6 


12 


60 


5.6 


69 


7 


9 


70 


4 


80 


8 


7 


81 


3 


92 


9 


6 


91 


2.2 


104 


10 


5 


100 


1.75 


115 



massive than for planetesimal radius of 100 km. These results 
show that when including the gas drag of the envelope, the re- 
duction factor in the formation time also depends on the surface 
density of the disc. Note that the fractional decrease of the whole 
formation time is greater than the fractional increase of the core 
mass. 

In summary, this section focused on the study of the forma- 
tion of a giant planet with a full evolutionary code, assuming 
an oligarchic growth regime for the core. We first studied the 
enhancement of the effective capture radius due to the presence 
of the gaseous envelope. The atmospheric gas drag enlarges the 
capture cross-section, causing a reduction in the formation time 
but also a significant increment in the mass of the core. We also 
studied the impact of the size of accreted planetesimals in our 
model. We found that the whole formation process is acceler- 
ated when considering smaller planetesimals and that the reduc- 
tion factor depends on the surface density of the feeding zone. 
However, no substantial change in the final mass of the core was 
identified. 

4. Discussion and summary 

In this paper we examined the formation of a Jupiter-like object. 
Motivated by the work of Thommes et al. (120031 1, we selected an 
oligarchic growth regime for the core as the planetesimal accre- 
tion rate. The numerical simulations were made with the code 
presented in a previous paper by Benvenuto & Brunini (120051 1, 
which was updated in this study (see Sect. |2). This code was 
developed for a self-consistent calculation of giant planet for- 
mation, so the growth of the core is coupled to the growth of the 
envelope. In this sense, the evolution of the envelope's density 
profile is a natural outcome of the code and a relevant quan- 
tity for an accurate estimation of the enhancement in the proto- 
planet capture cross-section of planetesimals (as the drag force 
depends on p s , see Eq. (fT51)), which has a direct impact on the 
growth rate of the core. 

When including a time-dependent planetesimal accretion 
rate in these kinds of calculations, most authors prescribe a rapid 
one, that in general guarantees the formation of a massive core 
in a few hundred thousand years. However, N-body simulations 
show that when the embryo is about the mass of the Moon, 
or even smaller (Thommes et al. 2003), the growth regime en- 
ters the protoplanet-dominated stage. For this reason, we as- 
sume from the very beginning of our calculations an oligarchic 
growth regime, which seems to be a more realistic approach 
to the problem. This accretion rate is considerably slower than 
the usually-adopted accretion model of Greenzweig & Lissauer 



(1992). Adopting the oligarchic growth model has, as a main 
consequence, an increase in the whole formation timescale when 
compared to previous calculations which adopt a rapid growth 
for the core (see Sect. 13. It . In this sense, it is worth mentioning 
that an accurate giant planet formation model with an oligarchic 
growth for the core should also include planetesimal migration 
due to the nebular gas drag since the assumed lifetime of the 
nebula is of the same order of magnitude of the formation pro- 
cess itself (Thommes et al. [2003 ). In the calculations presented 
here this effect was not included, as this article aims to analyse 
in a first approximation the oligarchic growth regime for a giant 
planet in a self-consistent calculation. 

To the authors' knowledge, the oligarchic growth of pro- 
toplanets has generally been studied in the absence of an at- 
mosphere and no full evolutionary calculations of gaseous gi- 
ant planets adopting this growth regime for the core have been 
performed in the past. As giant planets have large gaseous en- 
velopes, the gas drag effect of the growing atmosphere should 
play an important role in the resulting accretion cross-section of 
the protoplanet. In order to study the relevance of this effect, we 
made a set of simulations for several disc densities to estimate 
the consequences of including the gas drag in the calculation 
of the effective radius and we compared them to those where 
the effective radius is calculated in the absence of gas drag. 
According to our simulations, the main effects of the enlarge- 
ment of the effective radius when considering the gas drag are a 
reduction of about 30% to 55% in the whole formation timescale 
(depending on the surface density of the disc and on the plan- 
etesimal mass, being stronger for smaller planetesimals) and an 
increase in the mass of the core by about a factor of 2 for the 
several protoplanetary disc considered here. As we mentioned 
before, no planetesimal migration or other protoplanet-disc in- 
teractions were considered, so our results should be analysed in 
the context of this simple model. However, from the compari- 
son of simulations with and without atmospheric gas drag, we 
conclude that including the gas drag in the calculation of the ef- 
fective cross-section of the protoplanet is absolutely relevant to 
the estimation of the formation timescale. The fact that, in asso- 
ciation with higher surface densities for the disc, we obtain the 
desired shorter timescales but, on the other hand, super-massive 
cores should not be a cause for concern. Although the mass of 
the core is closely related to the formation timescale (the more 
massive the core, the shorter the timescale), Thommes et al. 
(2003) showed in their calculations of the formation of proto- 
planetary cores, that planetesimal migration strongly affects the 
planetesimal population of the feeding zone, being very effective 
in depleting it. Including this effect reduces considerably the fi- 
nal mass of the core, and protoplanetary discs as massive as 10 
MMSN would probably offer us good fits to the core mass and 
the formation timescale simultaneously. 

In this paper, we also explored the effect on giant planet for- 
mation of planetesimal size variations. For the same set of pro- 
toplanetary discs (from 6 to 10 MMSN), we made simulations 
for a swarm of planetesimals having a fixed radius of 10 and 100 
km. We find that the formation timescale is strongly dependent 
on planetesimal size: for planetesimal radius of 10 km, the pro- 
cess occurs a factor of 2 to 3 times faster than for the case of 
larger planetesimals. A reduction in the timescale was expected, 
since smaller planetesimals have lower relative velocities which 
favour accretion and they are much more affected by gas drag. 
However, the mass of the core is not very increased. 

In a recent paper, Inaba & Ikoma (2003) also explored the 
effects of gas drag on the effective radius of a protoplanet for 
capturing planetesimals of different sizes. They found that gas 
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drag largely increases the effective capture radius. While our re- 
sults are qualitatively similar to those of Inaba & Ikoma (2003), 
we note that we employed a full evolutionary code. In particu- 
lar, it is worth discussing the results Inaba & Ikoma (2003) pre- 
sented in their Fig. 6. They estimated the time spent in form- 
ing a 10 M ffi core as a function of planetesimal size, finding 
that gas drag drastically decreases the formation timescale: for 
planetesimals of, e.g., 10 km (100 km) the timescale is reduced 
by a factor of ~ 6 (~ 4). These factors would be very helpful 
in alleviating the timescale problem of the whole core growth 
mechanism. However, the results presented in this paper indicate 
that, while there exists a reduction in the time spent in form- 
ing a core of 10 M e , it is much more modest. We find that it 
is only slightly dependent on the protoplanetary surface density 
and moderately dependent on the planetesimal size. For plan- 
etesimals of 10 km the average reduction factor in the formation 
timescale of a 10 M ffi core is ~ 1.25 and for planetesimals of 
100 km it is ~ 1.5. It is also worth mentioning that the timescale 
involved in forming a 10 M e core is not really representative of 
the complete giant planet formation timescale as shown in this 
and previous papers (e.g., Pollack et al. 1 19961 . 

The next step for our future calculations is the inclusion of 
planetesimal migration due to nebular gas drag and planetesi- 
mal ejection out of the feeding zone. This will provide us with a 
more accurate quantitative idea of the final core masses and for- 
mation timescales in the frame of an oligarchic growth regime 
for the core. Also, a size distribution of the accreted planetes- 
imals should be considered for a better estimation of the for- 
mation timescales. Or, at least, if a unique radius is adopted for 
all planetesimals in the swarm, several simulations varying the 
radius should be made in order to bracket the formation time. 
Finally, we expect to update the opacity tables since opacity 
plays a fundamental role in the formation timescale (Podolak 
2003; Hubickyj et al. |20051 l. It would also be interesting to test 
with our model, the response of the whole planetary formation 
process to grain opacity reduction, which will probably shorten 
the timescales involved, as has been already shown by, e.g., 
Pollack et al. ( fl996l l and Hubickyj et al. (|2005I 
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